#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2018/1/31 16:56
# @Author  : Aries
# @Site    : 
# @File    : pms_spectra.py
# @Software: PyCharm Community Edition

from Py6S import *

def calc6SAtmosCorr(aeroProfile, atmosProfile, groundRefl, datestr, sensor_dic, surfaceAltitude, radi, useBRDF, aot550):
	# calculate 6s coeffs for each LUT combination
    reflarr = [-99, -99, -99, -99]  # 6s coeffs for 4 bands(red green blue n-vir)
    # s = SixS()
    s.atmos_profile = atmosProfile
    s.aero_profile = aeroProfile
    s.ground_reflectance = groundRefl

    s.geometry = Geometry.User()
    s.geometry.solar_z = 90 - sensor_dic['solar_zenith']
    s.geometry.solar_a = sensor_dic['solar_azimuth']
    s.geometry.view_z = 90 - sensor_dic['view_zenith']
    s.geometry.view_a = sensor_dic['view_azimuth']
    img_date = time.strptime(datestr, '%Y-%m-%d')
    day = img_date[2]  # datetime converts
    month = img_date[1]  # datetime converts
    s.geometry.day = day
    s.geometry.month = month
    s.geometry.gmt_decimal_hour = float(6) + float(01) / 60.0
    # s.geometry.from_time_and_location(31.6, 74.5, "14:24:59", 88.1613, 80.4918)

    # 气溶胶厚度使用能见度模式，默认为40 km，如有实测资料，也可使用550 nm 处的气溶胶厚度值来代替
    s.aot550 = aot550

    s.altitudes = Altitudes()
    s.altitudes.set_sensor_satellite_level()
    # s.altitudes.set_sensor_custom_altitude(645)
    # s.altitudes.set_target_sea_level()
    s.altitudes.set_target_custom_altitude(surfaceAltitude)

    if useBRDF:
        s.atmos_corr = AtmosCorr.AtmosCorrBRDFFromReflectance(radi)
    else:
        s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromRadiance(radi)

        # Band 1
    s.wavelength = returnGF1SpectraFilter(1)
    s.run()
    # print('b1 atmos_corrected_reflectance_brdf:', s.outputs.atmos_corrected_reflectance_brdf)
    # sixsCoeffs[0,0] = float(s.outputs.values['coef_xa'])
    # sixsCoeffs[0,1] = float(s.outputs.values['coef_xb'])
    # sixsCoeffs[0,2] = float(s.outputs.values['coef_xc'])
    reflarr[0] = s.outputs.atmos_corrected_reflectance_brdf

    # Band 2
    s.wavelength = returnGF1SpectraFilter(2)
    s.run()
    # print('b2 atmos_corrected_reflectance_brdf:', s.outputs.atmos_corrected_reflectance_brdf)
    # sixsCoeffs[1,0] = float(s.outputs.values['coef_xa'])
    # sixsCoeffs[1,1] = float(s.outputs.values['coef_xb'])
    # sixsCoeffs[1,2] = float(s.outputs.values['coef_xc'])
    reflarr[1] = s.outputs.atmos_corrected_reflectance_brdf

    # Band 3
    s.wavelength = returnGF1SpectraFilter(3)
    s.run()
    # print('b3 atmos_corrected_reflectance_brdf:', s.outputs.atmos_corrected_reflectance_brdf)
    # sixsCoeffs[2,0] = float(s.outputs.values['coef_xa'])
    # sixsCoeffs[2,1] = float(s.outputs.values['coef_xb'])
    # sixsCoeffs[2,2] = float(s.outputs.values['coef_xc'])
    reflarr[2] = s.outputs.atmos_corrected_reflectance_brdf

    # Band 4
    s.wavelength = returnGF1SpectraFilter(4)
    s.run()
    # print('b4 atmos_corrected_reflectance_brdf:', s.outputs.atmos_corrected_reflectance_brdf)
    # sixsCoeffs[3,0] = float(s.outputs.values['coef_xa'])
    # sixsCoeffs[3,1] = float(s.outputs.values['coef_xb'])
    # sixsCoeffs[3,2] = float(s.outputs.values['coef_xc'])
    reflarr[3] = s.outputs.atmos_corrected_reflectance_brdf
    return reflarr

def returnGF1SpectraFilter(band):
    """
	# read GF-1 pms&pan spectra data from file
	# parms: band type constants
	# return: spectra numpy array
	:param band:
	:return:
	"""
    # band type for choices
    # band = 'b0' for pan
    # band = 'b1' for B1
    # band = 'b2' for B2
    # band = 'b3' for B3
    # band = 'b4' for B4
    loc = 'I:\\scripts\\atmos_corr\\gf1_spectra.csv'  # location of spectra library in file system
    if band == 1:
        filterData = np.empty([30])
        i = 50
        j = 0
        while True:
            filterData[j] = npdata[i, 2]
            filterData[j + 1] = (npdata[i + 2, 2] + npdata[i + 3, 2]) / 2
            # print npdata[i,2]
            i = i + 5
            j = j + 2
            if npdata[i, 0] > 0.520:
                break
        # print 'b1:', filterData
        return Wavelength(0.4500, 0.5225, filterData)
    elif band == 2:
        filterData = np.empty([30])
        i = 120
        j = 0
        while True:
            filterData[j] = npdata[i, 3]
            filterData[j + 1] = (npdata[i + 2, 3] + npdata[i + 3, 3]) / 2
            # print npdata[i,3]
            i = i + 5
            j = j + 2
            if npdata[i, 0] > 0.590:
                break
        # print 'b2:', filterData
        return Wavelength(0.5200, 0.5925, filterData)
    elif band == 3:
        filterData = np.empty([26])
        i = 230
        j = 0
        while True:
            filterData[j] = npdata[i, 4]
            filterData[j + 1] = (npdata[i + 2, 4] + npdata[i + 3, 4]) / 2
            # print npdata[i,4]
            i = i + 5
            j = j + 2
            if npdata[i, 0] > 0.690:
                break
        # print 'b3:', filterData
        return Wavelength(0.630, 0.6925, filterData)
    elif band == 4:
        filterData = np.empty([50])
        i = 370
        j = 0
        while True:
            filterData[j] = npdata[i, 5]
            filterData[j + 1] = (npdata[i + 2, 5] + npdata[i + 3, 5]) / 2
            # print npdata[i,5]
            i = i + 5
            j = j + 2
            if npdata[i, 0] > 0.890:
                break
        # print 'b4:', filterData
        return Wavelength(0.770, 0.8925, filterData)
    else:  # pan band
        filterData = np.empty([182])
        i = 50
        j = 0
        while True:
            filterData[j] = npdata[i, 1]
            filterData[j + 1] = (npdata[i + 2, 1] + npdata[i + 3, 1]) / 2
            # print npdata[i,1]
            i = i + 5
            j = j + 2
            if npdata[i, 0] > 0.900:
                break
        # print 'pan:', filterData
        return Wavelength(0.450, 0.9025, filterData)

